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Abstract 

The  separation  of  variables  method  is  extended  to  obtain  concentration  profiles  in  a  particle  electrode  under  galvanostatic  boundary 
conditions.  The  method  is  also  used  to  find  exact  analytical  solutions  for  composite  slab  and  spherical  electrodes.  Finally,  the  method  is 
used  to  obtain  a  solution  for  a  lithium/polymer  cell  model  that  was  presented  previously  by  Doyle  and  Newman.  ©  2001  Elsevier  Science 
B.V.  All  rights  reserved. 
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1.  Introduction 

The  purpose  of  this  paper  is  to  present  an  extension  of  the 
separation  of  variables  method  for  solving  the  model  equa¬ 
tions  that  govern  concentration  distributions  in  solid  elec¬ 
trodes  operating  under  galvanostatic  conditions  (at  both  the 
ends).  The  method  presented  here  is  new  and  is  useful 
because  it  reduces  significantly  the  work  required  to  obtain 
an  analytical  solution  to  the  general  class  of  model  equations 
presented  here.  We  illustrate  our  method  for  a  thin  film 
electrode,  a  spherical  electrode  particle  and  composite 
electrodes. 


2.  Thin  film  electrode 


Consider  the  unsteady  state  diffusion  in  a  thin  film 
electrode  with  zero  initial  concentration.  The  governing 
equation  for  the  concentration  in  dimensionless  form  is 


dc_c?c 
dt  dx2 


(1) 


condition  c(x,  0)  =  0  and  boundary  condi  - 
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and 

g(U)=<5  <3) 

where  3  is  the  dimensionless  current  density. 

The  analytical  solution  for  this  boundary  value  problem 
(BVP)  is  given  by  [1] 

c  =(5|t+^(3x2-  1)— :  exp(— «27t2f)  cos(nnx) 

(4) 

This  solution  can  be  obtained  by  applying  the  Laplace 
transformation  technique.  Unfortunately,  inverting  back  to 
time  is  very  difficult  and  time  consuming  as  shown  in 
Appendix  A. 

Fortunately,  this  BVP  can  be  solved  easily  by  our  exten¬ 
sion  of  the  method  of  separation  of  variables  by  using  the 
following  variable  transformation: 

c(x,t)  =  u(x,t)+w(x)+v(t)  (5) 

The  addition  of  the  term  v(t)  in  Eq.  (5)  is  important  and  is  not 
presented  elsewhere  to  our  knowledge.  Here  w(x)  satisfies 
the  inhomogeneous  boundary  conditions 


0 


(6) 
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Nomenclature  1 

A 

separation  constant  (Eq.  (15)) 

An 

coefficient  in  the  infinite  series 

B 

constant 

C 

dimensionless  concentration  (Sections  2  and 

3  and  Appendix  A)  concentration  (mol/cm3, 
Sections  4  and  5) 

Cl,  c2 

concentration  (mol/cm3) 

Co 

initial  concentration  (mol/cm3) 

Ci,  C2 

dimensionless  concentration 

Di,D2 

diffusion  coefficient  (cm2/s) 

F 

Faraday’s  constant  (96,487  C/g  equivalent) 

i 

applied  current  density  (A/cm2) 

J 

dimensionless  current  density 

L 

length  (cm) 

r 

radial  distance  (cm) 

R 

radius  of  the  particle  (cm) 

s 

Laplace  transform  variable 

1 

dimensionless  time  (s,  Sections  2  and  3  and 
Appendix  A)  (Sections  4—6) 

X 

distance  (cm) 

X 

dimensionless  distance 

y 

dimensionless  distance 

Greek  letters 

d 

dimensionless  current  density 

4> 

dimensionless  potential 

A 

independent  variable  (only  in  Section  4) 

An 

eigenvalue 

Oi,  e2 

dimensionless  concentration 

T 

dimensionless  time 

The  transformation  given  by  Eq.  (5)  changes  Eq.  (1)  to 


du  dv  _  <92m  d2w 
dt  dr  dx2  +  dx2 
Separating  the  variables,  we  get 
du  _  d2u 

~di~W 


(12) 


(13) 


dv 


(14) 


dr  dx2 

Since  v(f)  is  a  function  of  t  only  and  vv(x)  is  a  function  of 
only,  we  require  that 


dv  _  d2w  _  ^ 


(15) 


where  A  is  a  constant.  Typically,  when  A  is  0,  Eq.  (5)  can  be 
replaced  by  the  usual  variable  transformation  [2], 

c(x,  r)  =  m(x,  r)  +  w(x)  (16) 

For  this  case,  A  is  nonzero  and  is  determined  by  the 
boundary  conditions.  The  second  half  of  Eq.  (15)  can  be 
solved  with  the  boundary  conditions  to  give 


(17) 


and 


w(x)  =1<5x2+S  (18) 

where  B  is  an  arbitrary  constant.  The  left-hand  side  of 
Eq.  (15)  can  be  solved  with  the  initial  condition  for  v(r) 
(Eq.  (10))  to  give 


and 


v(r)  =  dt 


(19) 


7(1) 


(7) 


The  variable  u(x,  t )  satisfies  the  homogeneous  boundary 
conditions 


du 

dx 


(0,0  =  0 


and 


(8) 


The  variable  v(r)  satisfies  the  initial  condition 
v(0)  —  0 


(9) 


(10) 


and  all  three  variables  together  also  satisfy  the  initial  con¬ 
dition 


Hence,  the  solution  is  given  by 

c(x,  r)  =  u(x,  t)  +  w(x)  +  v(r)  =  u(x,  r)  + 1  dx1  +  dt  +  B 

(20) 

Now  m(x,  r)  is  obtained  by  solving  Eq.  (13)  with  the  homo¬ 
geneous  boundary  conditions  (Eqs.  (8)  and  (9))  to  give 

m(x,  r)  =  y>,  exp(— «27t2r)  cos(rarx)  (21) 

where  A„  (n  -  1.2,...)  are  constants.  Hence,  the  final 
solution  is  given  by 

c(x,  t)  =  ^x2  +  dt  +  B  +  yAn  exp(—n2n2t)  cos(nnx) 

(22) 

The  constants  A„  and  B  are  obtained  by  imposing  the  initial 
condition 


c(x,  0)  =  0  =  -x2  +  B  +  YAn  cos(nnx) 


((x,  0)  +  w(x)  +  v(0)  =  0 


(ID 


(23) 
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Eq.  (23)  is  of  Sturm-Liouville  type.  A„  can  be  obtained  by 
multiplying  both  sides  of  Eq.  (23)  by  cos(nnx)  and  integrat¬ 
ing  from  0  to  1  to  get 


The  constant  B  can  be  obtained  from  Eq.  (23)  by  specifying  a 
value  for  x  (and  picking  10  or  more  terms  in  the  expansion) 
and  solving  for  B.  Alternatively,  B  can  be  obtained  by 
integrating  both  sides  of  Eq.  (23)  from  0  to  1 

B=-\d  (25) 

Substituting  Eqs.  (24)  and  (25)  into  Eq.  (22)  yields  Eq.  (4). 

Mathews  and  Walker  [2]  solved  the  same  problem  by 
using  the  transformation  given  by  Eq.  (16).  One  cannot  solve 
explicitly  for  w(x)  by  using  Eq.  (16)  subject  to  the  boundary 
conditions  given  by  Eqs.  (6)  and  (7).  They  assumed  arbi¬ 
trarily  a  parabolic  profile  for  w(x )  satisfying  these  two 
boundary  conditions  (Eqs.  (6)  and  (7)),  thereby  make  the 
governing  equation  for  u(x,  t )  inhomogenous.  Consequently, 
their  method  leads  to  a  more  difficult  problem  to  solve  and  is 
not  general. 

Clearly  our  variable  transformation  Eq.  (5)  is  a  more 
direct  method  than  that  presented  by  Mathews  and  Walker 
(we  do  not  have  to  assume  that  w(.x)  is  parabolic)  and  our 
method  is  easier  to  apply  than  the  classical  Laplace  trans¬ 
form  technique  given  in  Appendix  A. 

A  similar  problem  is  solved  by  Bird  et  al.  [3]  (pages  295- 
296  and  362)  for  laminar  tube  flow  with  constant  heat  flux  at 
wall.  However,  they  assumed  the  form  of  the  solution  (see 
Eq.  (9.8.23)  of  [3]).  Also,  they  introduced  another  boundary 
condition  (Eq.  (9.8.25)  of  [3])  by  heat  balance.  Our  method 
does  not  require  assuming  the  form  of  the  solution  as  in 
references  [2,3],  Also,  our  method  does  not  need  the  addi¬ 
tional  boundary  condition  introduced  by  [3]. 


3.  Spherical  electrode  particle 


The  utility  of  our  method  is  even  greater  for  the  case  of  a 
spherical  electrode  particle  with  unsteady  state  diffusion  and 
zero  initial  concentration.  The  governing  equation  for  the 
concentration  in  dimensionless  form  is  given  by 


dc  1  <9  /  2*5c\ 
~dt=^dr\dr) 


(26) 


with  the  initial  condition  c(r,  0)  =  0  and  boundary  condi¬ 
tions 


|«M)  =  0  (27) 

and 


^<1,0  =  <5  (28) 

where  again  (5  is  the  dimensionless  current  density. 


The  analytical  solution  for  this  BVP  is  given  by  [4] 


(29) 


where  ln  {n  =  1 , 2, . . .)  are  the  positive  roots  of 
tan(2„)  =  (30) 

This  solution  can  be  obtained  by  applying  the  Laplace 
transformation  in  the  time  variable,  as  shown  in  Appendix 
A.  Unfortunately,  using  the  Laplace  transform  technique  for 
this  problem  is  even  more  difficult  and  time  consuming. 
Fortunately,  our  method  works  well  and  is  easy  to  apply.  As 
before,  the  following  variable  transformation  is  used: 

c(r,  t)  =  u(r ,  t)  +  w(r)  +  v(t)  (31) 

where  w(r)  satisfies  the  inhomogeneous  boundary  condi¬ 
tions 
dvr 

-(0)=0  (32) 

and 


The  variable  u(r,  t)  satisfies  the  homogeneous  boundary 
conditions 


dr 

and 


dii 


0 


(35) 


The  variable  v(t)  satisfies  the  initial  condition 
v(0)  =  0  (36) 

and  three  variables  together  also  satisfy  the  initial  condition 
u(r,0)  +  w(r )  +  v(0)  =  0  (37) 

Applying  the  same  procedure  as  before  the  variables  are 
solved  as 


w(r)  =  t5r>- 
v(t)  =  3<5t 
and 


(38) 

(39) 


u(r,  t)  =  sin(A„r)  exp (— 


(40) 


where  A„  (n  -  1.2,...)  are  constants.  Hence,  the  final 
solution  is  given  by 

c(r,  t)  =  *^r1  +  5t  +  B  +  - y^A„  sin(2„r)  exp(— )?nt)  (41) 
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The  constants  An  and  B  are  obtained  by  imposing  the  initial 
condition 


c(r,0)  =  0  =  ^r2 


-y>„  sin(2„r) 


(42) 


Eq.  (42)  is  of  Sturm-Liouville  type.  A„  can  be  obtained  by 
multiplying  both  sides  of  Eq.  (42)  by  r  sin(A„r)  and  inte¬ 
grating  from  0  to  1  to  get 


An  = 


2 

sin(A„) 


(43) 


The  constant  B  is  obtained  by  multiplying  both  sides  of 
Eq.  (42)  by  r2  and  integrating  from  0  to  1 

B=~^5  (44) 


4.  Composite  electrodes 

We  consider  two  types  of  composite  electrodes  in  this 
section.  One  is  a  thin  film  composite  electrode  with  three 
layers.  The  inner  layer  is  made  of  material  different  from  the 
outer  two  layers.  An  electrochemical  reaction  is  assumed  to 
be  occurring  on  the  surfaces  of  the  outer  layers  so  that  the 
centerline  of  the  inner  layer  can  be  considered  to  be  a  plane 
of  symmetry.  The  second  composite  electrode  (Section  5)  is 
a  spherical  electrode  with  an  inner  core  and  an  outer  layer 
where  an  electrochemical  reaction  is  assumed  to  occur  at  the 
surface.  These  composite  electrodes  can  be  used  to  model 
intercalation  electrodes  found  in  the  nickel/metal  hydride 
and  lithium  ion  batteries  where  charge  is  stored  in  a  solid 
phase  in  the  form  of  charged  species  (e.g.  Li+).  Unfortu¬ 
nately,  no  analytical  solution  exists  in  the  literature  (that  we 
could  find)  for  a  composite  electrode  with  flux  conditions. 

Fortunately,  our  method  can  be  used  to  derive  an  analy¬ 
tical  solution  for  such  cases.  Consider  a  rectangular  slab 
with  region  1  (0  <  x  <  l\)  of  one  medium  with  diffusion 
coefficient  Di  and  region  2  (/,  <  x  <  L)  of  another  medium 
with  diffusion  coefficient  D2,  respectively.  The  governing 
equations  for  the  concentration  of  diffusing  species  in 
regions  1  and  2  (ci  and  c2)  are 


»<*<'■ 

(45) 

siF 

b 

gj|  ^ 

A 

A 

(46) 

Let  the  initial  concentration  be  c0,  i.e. 

ci(x,0)  =  c2(x,  0)  =  c0,  for  all x 

(47) 

Symmetry  at  the  center  (x  =  0)  gives  the  boundary  condition 

(0,  t)  =0,  for  t  >  0 
dx 

(48) 

Under  galvanostatic  discharge  conditions,  applied  current 
gives  the  boundary  condition  at  the  surface 

d c 2 


fort  >  0 


(49) 


where  i  is  the  applied  current  density  (a  positive  quantity)  at 
the  surface  of  the  electrode,  F  the  Faraday’s  constant  and  n 
the  number  of  electrons  transferred  in  the  electrochemical 
reaction.  The  negative  sign  in  Eq.  (49)  is  used  because 
charge  is  being  removed  from  the  battery.  Assuming  no 
charge  transfer  resistance  at  the  interface  (x  =  li),  the 
boundary  conditions  at  the  interface  are 


>  0 


ci(/i,f)  =  c2(h,t), 

;>o 


dx  y 


(50) 

(51) 


Next,  we  introduce  the  following  dimensionless  variables: 

Ci  —  —  —  i-  C2  =  «- 1,  x  =  *  t  = 

Co  Co  L  Lr 

(52) 

The  governing  equations  in  dimensionless  form  are 

0  <  X  <  a  (53) 


dCi  1  (fC\ 
dx  ~  p2  dX2  ’ 


dC2_ 

dx 

where 

h 


d2C2 

dx2  ’ 


and 

F  =  %  (56) 

L>\ 

The  initial  and  boundary  conditions  are  transformed  to 
Ci  (X,  0)  =  C2  (X,  0)  =  0,  for  all  X  (57) 


dX  y 

Ci  (a,  t)  =  C2(a,  x)  =0,  x  >  0 
dCi.  2  dC2 

«(“■')  =  '*  t>0 


(58) 

(59) 

(60) 
(61) 


<5  = 


iL 

nFD2c0 


(62) 


is  the  applied  dimensionless  current  density  at  the  surface  of 
the  electrode  (a  positive  number). 

Now  applying  our  variable  transformation  defined  by 
Eq.  (5)  for  both  Ci  and  C2  separately  and  solving  with 
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the  initial  (Eq.  (57))  and  boundary  conditions  (Eqs.  (58)  and 
(59)),  we  get 


Ci  =  ]- kf}2X 2  +  kx  +  a\  +  y ^Bn  cos(X„px)  exp(-^t): 


C2 


X-kX 2  -(6  +  k)X  +  kx  +  a2  +  cos(A„[l  -  X]) 

x  exp(— A2t)  (64) 


where  k,  a\  and  a2  are  constants  to  be  determined  by  the 
boundary  conditions  at  the  interface  ( X  =  a),  Xn 
(n=  1,2,3,...)  are  the  eigenvalues,  and  Bn  and  En  are 
the  constants  associated  with  the  eigenvalues.  Boundary 
conditions  Xn  at  the  interface  (Eqs.  (60)  and  (61))  yield 

k  =  -5  (65) 

a2  —  0i  =  |<5a2(l  —  /?2)  (66) 

B„  En 

cos(A„[l  -  a])  cos(A„/Ja)  A" 

where  A„  is  a  new  constant  associated  with  the  eigenvalue  Xn. 
Now  the  solutions  for  Ci  and  C2  are  given  by  the  simplified 
expressions 


Cl  =  -  ljdf}2X2  -  St  +  a\  +  y~/l„  cos(A„[l  -  a]) 

2  «=  l 

x  cos (Xnpx)  exp(— A2t)  (68) 

C2  =  —  ^c)X2  —  Sz  +  02  +  y>„  cos(A„[l  —  X]) 

x  cos(A„/?a)  exp(— X2z)  (69) 


where  Xn  are  the  positive  roots  of  the  transcendental  equation 
,»(!„[! -«])  +  ,a^  =  0  ,70) 


film  electrode,  the  constants  are  obtained  as  follows: 

«i  =  ^(s  +  5«3[1  -  (52]  -|«2[1  -  >?2])  (73) 

02  =  <5(§  +  5«3[l-£2])  (74) 

and 

A  = _ 28cfs{$\*) _ 

acos2(A„[l  —  a])  +  (1  —  a)  cos2(/Ll„a) 

Hence,  the  concentrations  of  the  diffusing  species  are  given 
by 

Cl  =  Co  |l  -  p2X 2  -  St  +  ai  +  ^A,,  cos(A„[l  -  a]) 

x  cos(X„[5X )  exp(— A2t)J 

c2  =  c0 1 1  —  ^<SX2  —  Sz  +  a2  +  y,A„  cos(A„[l  —  X]) 

x  cos(A„/?a)  exp(-A2ir)j  (77) 

where  a\,a2  and  A„  are  defined  by  Eqs.  (73)-(75)  and  the 
eigenvalues  X„  are  given  by  Eq.  (70).  When  /?  =  1,  both  ci 
and  c2  given  by  Eqs.  (76)  and  (77)  reduce  to 

c  =  c0  |l  -  <5t  -  ^  (3X2  exp(— n27r2i) 

x  cos(n7tX)j  (78) 

which  is  the  solution  for  diffusion  in  a  thin  film  electrode 
of  diffusion  coefficient  D2  ( D2  =  D\  =  D),  as  expec¬ 
ted.  Dimensionless  concentration  (ci/co[0  <  X  <  a]  or 
c2/co[a  <  X  <  1])  profiles  for  5  =  1,  a  =  0.5  and  P  =  2 
are  given  in  Fig.  1  for  several  values  of  t. 


5.  Composite  spherical  electrodes 


which  was  obtained  by  equating  the  flux  at  the  interface 
(Eq.  (61)).  Now  Ci  and  C2  separately  satisfy  the  initial 
condition  (Eq.  (57)) 

Ci (X,  0)  =  0  =  -  X- dp2X2  +ai  +  cos{Xn  [1  -  a]) 

x  cos (Xnpx)  (71) 

and 

C2(X,  0)  =  0  =  —  ^SX2  +  a2  +  yAn  cos (An [1  -  X]) 

x  cos(XJa)  (72) 

Eqs.  (71)  and  (72)  are  of  the  Sturm-Liouville  type.  Applying 
a  procedure  similar  to  the  one  presented  above  for  the  thin 


A  general  model  for  a  composite  spherical  electrode  is 
presented  here  by  using  our  approach.  The  governing  equa¬ 
tions  for  the  concentrations  of  the  diffusing  species  are 


dc\  Di  d  (  ,  dcA 

(79) 

~di~Wd?{r  ^ )’  0<r<n 

dc2  D2  d  (  ,  dc2\ 

Ot  =  r2  dr  \  dr)'  r'<r<R 

(80) 

with  the  initial  and  boundary  conditions 

ci(r,  0)  =  c2(r,  0)  =  Co,  for  all  r 

(81) 

■^(0,0  =  0,  t>o 

(82) 

(83) 
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ci (n ,  t)  =  c2(n ,  t) ,  t  >  0  (84) 

=  D2-g^(n,t),  t>  o  (85) 


Using  the  same  procedure  followed  for  a  composite  rectan¬ 
gular  slab  presented  above,  we  arrive  at  the  solution  for  the 
concentration  profiles 


Cl  =  CO  [l  -  \w2X2  -  3dz 

+  ai  +x^An 

x  [A„  cos(A„[l  —  a])  — 

sin(A„[l  —  a])]  sin 

(W 

x  exp(— A^r)| 

(86) 

ci  =  c0  |l  -  ^(3X2  -  3St  + 

1  ^ 

x  [A„cos(A„[l  —  X]  )  - 

■  sin(A„[l  —  X])]  sin(A„/?a) 

x  exp(-A^r)] 

(87) 

where 

5=  lR 

nFDiCo 

(88) 

X-r  z~Dlt 

X  R'  R 2 

(89) 

n 

R 

(90) 

P  DX 

(91) 

a]=S(l+la5[)-p2]-l 

^[l-p2]) 

(92) 
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and  the  boundary  conditions  (Eqs.  (106)-(109))  become  Thus,  the  dimensionless  concentration  profiles  are 


dC,  r  n 

ay 

(115) 

0i  =  1  +  J  +  er  +  cos  cos(A„y)  exp(— A^t) j  , 

dSl=0*y=l  +  r 

(116) 

Ci  =  C2  aty  =  1 

(117) 

+  fSl„cos^  '  gl/4  cos(A„)  exp(  A^t)J 

^i=^at5tai 
ay  ay 

(US) 

(126) 

with  the  initial  conditions 

Ci  =  C2  =  0  at  r  =  0  (119) 

This  problem  is  very  similar  to  the  composite  electrode 
problem  solved  in  Section  4.  We  apply  the  variable  trans¬ 
formation  (Eq.  (5))  for  both  Ci  and  C2 

Ci  =  Mi(y,r)  +  wi(y)  +  vi(t)  (120) 

and 

C2  =  u2(y,  t)  +  w2(y)  +  v2(r)  (121) 

By  applying  the  same  procedure  as  before  we  get 
Vi(r)  =  v2(r)  =  0,  Wi(y)=B  +  sry, 

Wl(y)  =  B  +  ,ry-±i(l  +  r)-±;^-y(l  +  r)) 

(122) 

Note  that  vtq  and  w2  satisfy  all  the  four  boundary  conditions 
given  by  Eqs.  (1 15)— (1 18).  The  transient  solutions  are 
obtained  as  before 


where  the  constants  B  and  A„  are  given  by  Eq.  (125).  Note 
that  Eq.  (126)  is  simpler  than  the  solution  reported  by  Doyle 
and  Newman  [5],  The  solution  given  by  Doyle  and  Newman 
involves  two  different  constants  ( Fn  and  G„)  for  both  9i  and 
92  separately.  Even  though  our  solution  looks  different  from 
that  of  Doyle  and  Newman,  both  are  equivalent.  The  con¬ 
centration  is  minimum  at  the  end  of  the  cathode  (y  =  1  +  r). 
The  time  for  this  concentration  to  reach  zero  is  the  dimen¬ 
sionless  depletion  time.  This  concentration  can  be  obtained 
by  substituting  y  =  1  +  r  into  Eq.  (114).  Doyle  and  Newman 
used  only  one  term  in  the  expansion  for  longer  times.  They 
found  a  short  time  solution  by  applying  the  Laplace  trans¬ 
form  technique.  It  should  be  noted  that  as  mentioned  by 
Doyle  and  Newman  one  term  in  the  series  in  Eq.  (14)  and  the 
short  time  solution  is  sufficient  for  estimating  the  depletion 
time.  However,  if  one  wants  to  follow  the  transient  behavior 
of  the  cell  sandwich,  one  term  is  not  sufficient.  We  have  used 
N  =  20  terms  and  predicted  the  concentration  profiles  inside 
the  cell  sandwich  for  various  values  of  time  and  plotted  them 
in  Fig.  3  for  the  given  value  of  applied  current,  J  =  — 0. 1 .  As 
shown  in  Fig.  3,  the  one  term  approximation  used  by  Doyle 
and  Newman  is  not  valid  for  predicting  the  transient  beha¬ 
vior  inside  the  cell  sandwich. 


mi  =  f>„cos(^)  cos(A„y)  exp(— a^t), 

m2  =  cos  ^4 — — ^  cos(A„)  exp(-A^i)  (123) 


where  A„  (n  =  1 , . . . ,  oo)  and  the  eigenvalues  are  given  by 
the  transcendental  equation 

tan(A„)  =  -e5/4  tan(e“1/4  A„r)  (124) 

By  applying  the  initial  condition  and  by  using  the  Sturm- 
Louiveille  theorem  (see  Section  4)  we  get 


B  = 
An= 


1  /fir  g2?.2\  1  Vzf3 

1  +sr  V  2  8  Y  )  3  1+  er  ’ 

2A„  cos (A„r/e1/4)£r  —  e5/4  sin(2„r/e1/4)cos(2„) 
A;)  (cos2  (A„r/ e1/4) +er  cos2  (A„) ) 


(125) 


c 


Fig.  3.  Dimensionless  transient  concentration  profiles  in  a  time  in  a 
lithium  ion  cell  sandwich. 
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Fig.  4.  Concentration  profiles  in  a  lithium  ion  cell  sandwich  as  a  function 
of  applied  current,  J  at  a  particular  time  (t  =  1). 

We  would  like  to  stress  again  that  we  obtained  a  simpler 
but  equivalent  solution  to  that  obtained  by  Doyle  and  New¬ 
man.  The  effect  of  applied  current,  J  can  be  seen  using  the 
same  set  of  coefficients  with  our  method  unlike  Doyle  and 
Newman.  This  is  true  because  their  coefficients  are  functions 
of  the  applied  current,  J  (Eqs.  (33)  and  (34),  [5])  which  means 
that  one  has  to  calculate  the  coefficients  for  every  value  of  J. 
At  a  particular  value  of  dimensionless  time  (t  =  1),  the 
concentration  profiles  inside  the  cell  sandwich  are  plotted 
in  Fig.  4  for  different  value  of  applied  current  density,  J. 


1.  Conclusion 

An  extension  of  the  method  of  separation  of  variables 
presented  here  and  given  by  Eq.  (5)  is  useful  for  solving 
BVPs  that  include  flux  boundary  conditions.  This  method 
yields  an  unambiguous,  straightforward  way  to  obtain  ana¬ 
lytical  solutions  for  problems  of  this  type.  The  utility  of  the 
method  is  demonstrated  for  two  classical  problems  (diffu¬ 
sion  in  a  slab  and  a  sphere).  Also  the  method  is  used  to  obtain 
analytical  solutions  for  diffusion  in  slab  and  spherical 
composite  electrodes.  The  transformation  proposed  helps 
in  obtaining  a  compact,  and  easy  to  use  solution  for  the 
lithium  ion  cell  sandwich  under  solution  phase  diffusion 
limitations.  The  method  applied  can  be  easily  extended  to 
cylindrical  coordinates.  The  method  appears  to  be  general 
and  should  be  useful  for  solving  other  similar  problems  with 
flux  boundary  conditions. 
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Appendix  A.  Thin  film  electrode 

The  Laplace  transform  (for  the  time  variable)  can  be 
applied  to  Eqs.  (l)-(3)  to  obtain 


d2c(s) 

dx2 


(A.l) 
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Substituting  5  =  0  and  simplifying  we  get 

m  =  s 


(A.  12) 
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Similarly  from  Eq.  (A. 8) 


(j)'{s)  = 


5  cosh(-v/sx)  sinh(-v/5;)  +  Xyfs  sint^y^x)  sinh(^/i)  —  -v/5:cosh(-v/ix)  cosh(-v/5) 

2  V^[sinh2(V^)] 


(A.  13) 


which  for  s  =  0  yields 


Eq.  (A.  19)  can  be  written  in  series  form 


0'(O)  =  (A.  14) 

After  applying  L’  Hospital’s  rule  three  times  we  get 

=  +  ^  (A- 15) 

(The  expressions  obtained  after  applying  L’  Hospital’s  rule 
are  too  big  to  be  included  here.)  The  complete  solution  can 
be  obtained  by  using  Heaviside  expansion  theorem  for  no 
repeated  roots  (Eq.  (8.3.22)  of  [6])  and  adding  the  result  to 
Eq.  (A.9);  the  solution  is  given  by 


_  P(s)  1  o^("+1)/2r2"+1/(2n  +  1)! 

1  J  Q(s)  +  l)(2n  -  1)! 

(A.20) 

In  Eq.  (A.20),  the  polynomial  Q(s )  has  a  greater  power  in  s 
than  P(s).  The  poles  are  s  =  0  (multiplicity  2)  and  the  roots 
of  tan(2„)  =  Xn(tt=  1, 2, 3, . . .),  where  s  =  — /2.  Using  the 
Heaviside  expansion  theorem  for  repeated  roots  [6] 


d>(s)  =  (s-a)mf(s)  = 


5s  sinh(\/yr) 
y/s  cosh  y/s  —  sinh(  y(s) 


(A.21) 


c(x, t)  =  St-^  +  ^x2  +  |>p(  -  l2nt)  (p(2„)  (A.  16) 

where  Xn  =  nn  and 


From  Eq.  (A.21)  we  have  (for  5  =  0) 


m 


(*)(Q)Q  _  o 

r(0  -  0)  0 


(A.22) 


cp(Xn)  =(s-  Xn)f(s )  = 


P(s  =  ~K) 

®f  =  ~%) 


Hence,  the  complete  solution  is 
x  exp(  —  n2n2t)  cos(nra;)J 


—25  cos(nnx) 

(A.  17) 


(A.  18) 


After  applying  L’  Hospital’s  rule  to  Eq.  (A.21)  two  times  and 
simplifying 

=  <j3cosh('/^)+W»inh(v%) 
cosh  y/s 

Substituting  5  =  0  and  simplifying  we  get 

0(0)  =  35  (A.24) 

Similarly  from  Eq.  (A.21)  we  get 


5  2 y/s  sinh(-v/5x)  cosh(y/5)  —  (5  +  2)  sinh(-v/5x)  sinh(y/5)  +  xs  cosh(^/5x)  cosh(\/5)  —  Xy/s  cosh(-v/5x)  sinh(-v/5) 

2  5COsh2(y/5)  —  2y/s cosh(-v/5)  sinh(y/5)  +  sinh2(^/5) 

_  (A.25) 


As  seen  for  this  case,  inversion  to  the  time  domain  is 
cumbersome  even  for  simple  eigenvalues  of  Xn  =  nn.  Eva¬ 
luation  of  the  residues  at  5  =  0  involved  using  L’  Hospital’s 
rule  thrice  even  for  this  simple  case. 


A.l.  Spherical  electrode  particle 

As  before,  applying  the  Laplace  transform  (in  the  time 
variable)  to  Eq.  (26)  and  solving  with  the  boundary  condi¬ 
tions  (27)  and  (28)  we  get 


1  5  sinhfi/yr)  P(s) 

r  s(y/s  coshyfs  —  sinh(-v/5))  Q(s) 


(A.  19) 


which  for  5  =  0  yields 

0'(O)  =  ?  (A.26) 

After  applying  L’  Hospital’s  rule  six  times  we  get 
0*(O)  =  -|:|^J  Sr2  (A.27) 

When  we  go  from  rectangular  to  spherical  coordinates  even 
for  finding  the  residue  at  the  origin,  we  need  to  L’  Hospital’s 
rule  twice  the  number  of  times.  Next,  by  using  the  Heaviside 
expansion  theorem  for  no  repeated  roots  (Eq.  (8.3.22)  of 
[6]),  the  solution  in  the  time  domain  is  given  by 

c(x,  t)  =  35t  -  ^  +  ^r2  +  |]/xp (-A2f)<p(2„)  (A.28) 
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where  Xn  is  given  by  tan(A„)  =  a„  and 
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Q'(s  =  -^) 


1  2 S  sin(V) 

r  sin(2„) 

(A.29) 


Hence,  the  complete  solution  is 


m 


(A.30) 
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